function [NewEqtn,NewEqtnF,NewEqtnS,IxNewNonlin] ...
    = myoptpolicy(This,PosLossEqtn,LossDisc,PosNnegName,PosNnegMult,Type)
% myoptpolicy  [Not a public function] Calculate equations for discretionary optimal policy.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

%                      | Original   | With optimal    | With optimal
%                      |            |                 | and nonneg
%----------------------+------------+-----------------+---------------------
%                      |            |                 | 
%                      | Eq1   Var1 | Eq1:=dMu1  Mu1  | Eq1           Mu1
%                      | Eq2   Var2 | Eq2:=dMu2  Mu2  | Eq2           Mu2
% PosLossEqtn          | Loss  Var3 | dVar1      Var1 | Min           MuNneg
%                      |       Var4 | dVar2      Var2 | dVar1         Var1
%          PosNnegName |       Var5 | dVar3      Var3 | dVar2=MuNneg  Var2
%                      |            | dVar4      Var4 | dVar3         Var3
%                      |            | dVar5      Var5 | dVar4         Var4
%                      |            |                 | dVar5         Var5

%--------------------------------------------------------------------------

isNneg = ~isempty(PosNnegName);
template = sydney();
flNameType = floor(This.nametype);
t0 = find(This.Shift==0);

eqtn = cell(size(This.eqtnF));
eqtn(:) = {''};
eqtn(flNameType==2) = This.eqtnF(flNameType==2);

% Replace x(:,n,t+k) with xN, xNpK, or xNmK, and &x(n) with Ln.
eqtn = sydney.myeqtn2symb(eqtn);
LossDisc = sydney.myeqtn2symb(LossDisc);

% First transition equation.
first = find(This.eqtntype==2,1);

% The loss function is always the last equation. After the loss function,
% there are empty place holders. The new equations will be put in place of
% the loss function and the placeholders.
NewEqtn = cell(size(eqtn));
NewEqtn(:) = {''};
IxNewNonlin = false(size(eqtn));
zd = sydney(LossDisc,{});

% The lagrangian is
%
%     Mu_Eq1*eqtn1 + Mu_Eq2*eqtn2 + ... + lossfunc.
%
% The new k-th (k runs only for the original transition variables) equation
% is the derivative of the lagrangian wrt to the k-th variable, and is
% given by
%
%     Mu_Eq1*diff(eqtn1,namek) + Mu_Eq2*diff(eqtn2,namek) + ...
%     + diff(lossfunc,namek) = 0.
%
% We loop over the equations, and build gradually the new equations from
% the individual derivatives.

for eq = first : PosLossEqtn
    
    % Get the list of all variables and shocks (current dates, lags, leads) in
    % this equation.
    [tmOcc,nmOcc] = myfindoccur(This,eq,'variables_shocks');
    tmOcc = tmOcc(:).';
    nmOcc = nmOcc(:).';

    % Find the index of transition variables.
    inx = flNameType(nmOcc)==2;
    if strcmpi(Type,'consistent') || strcmpi(Type,'discretion')
        % This is a consistent (discretionary) policy. We only
        % differentiate wrt to current dates or lags of transition
        % variables. Remove leads from the list of variables we will
        % differentiate wrt.
        inx = inx & tmOcc <= t0;
    end
    
    tmOcc = tmOcc(inx);
    nmOcc = nmOcc(inx);
    nOcc = length(tmOcc);
    
    % Write a cellstr with the symbolic names of variables wrt which we will
    % differentiate.
    unknown = cell(1,nOcc);
    for j = 1 : nOcc
        if tmOcc(j)==t0
            % Time index==0: replace x(1,23,t) with x23.
            unknown{j} = sprintf('x%g',nmOcc(j));
        elseif tmOcc(j)<t0
            % Time index<0: replace x(1,23,t-1) with x23m1.
            unknown{j} = sprintf('x%gm%g', ...
                nmOcc(j),round(t0-tmOcc(j)));
        elseif tmOcc(j)>t0
            % Time index>0: replace x(1,23,t+1) with x23p1.
            unknown{j} = sprintf('x%gp%g', ...
                nmOcc(j),round(tmOcc(j)-t0));
        end
    end
    
    z = sydney(eqtn{eq},unknown);
    
    % Differentiate this equation wrt to all variables and return a cell array
    % of separate sydney objects for each derivative.
    diffz = derv(z,'separate',unknown);
  
    for j = 1 : nOcc

        sh = tmOcc(j) - t0;
        newEq = nmOcc(j);
        
        % Multiply derivatives wrt lags and leads by the discount factor.
        if sh==0
            % Do nothing.
        elseif sh==-1
            diffz{j} = zd * diffz{j};
        elseif sh==1
            diffz{j} = diffz{j} / zd;
        else
            diffz{j} = power(zd,-sh) * diffz{j};
        end
        
        % If this is not the loss function, multiply the derivative by the
        % respective multiplier. The appropriate lag or lead of the multiplier will
        % be introduced together with other variables.
        if eq<PosLossEqtn
            mult = template;
            mult.args = sprintf('x%g',eq);
            diffz{j} = diffz{j}*mult;
        end

        dEqtn = char(reduce(diffz{j}));
        
        % Shift lags and leads of variables and multipliers (but not parameters) in
        % this derivative by -sh if sh~=0.
        if sh~=0
            dEqtn = sydney.myshift(dEqtn,-sh,flNameType<=2);
        end
        
        sign = '+';
        if isempty(NewEqtn{newEq}) ...
                || strncmp(NewEqtn{newEq},'-',1) ...
                || strncmp(NewEqtn{newEq},'+',1)
            sign = '';
        end
        NewEqtn{newEq} = [dEqtn,sign,NewEqtn{newEq}];
        if ~This.IsLinear
            % Earmark the derivative for non-linear simulation if at least one equation
            % in it is nonlinear and the derivative is nonzero. The derivative of the
            % loss function is supposed to be treated as nonlinear if the loss function
            % itself has been introduced by min#() and not min().
            isNonlin = This.IxNonlin(eq) && ~isequal(dEqtn,'0');
            IxNewNonlin(newEq) = IxNewNonlin(newEq) || isNonlin;
        end
        
    end
end

% Add nonnegativity multiplier to RHS of derivative wrt nonegative
% equation, and create complementary slack condition, which is always
% placed where loss function was.
if isNneg
    NewEqtn{PosNnegName} = [NewEqtn{PosNnegName},sprintf('-x%g',PosNnegMult)];
    NewEqtn{PosLossEqtn} = sprintf('min(x%g,x%g)',PosNnegName,PosNnegMult);
    IxNewNonlin(PosLossEqtn) = true;
end

NewEqtnF = sydney.mysymb2eqtn(NewEqtn,'full');
NewEqtnS = sydney.mysymb2eqtn(NewEqtn,'sstate');
NewEqtn = sydney.mysymb2eqtn(NewEqtn,'human',This.name);

% Add semicolons at the end of each new equation.
pos = PosLossEqtn : length(NewEqtn);
NewEqtn(pos) = strcat(NewEqtn(pos),'=0;');
NewEqtnF(pos) = strcat(NewEqtnF(pos),';');
NewEqtnS(pos) = strcat(NewEqtnS(pos),';');

% Replace = with #= in nonlinear human equations.
if ~This.IsLinear
    NewEqtn(IxNewNonlin) = strrep(NewEqtn(IxNewNonlin),'=0;','=#0;');
end
    
end
